Density wave instabilities of tilted fermionic dipoles 
in a multilayer geometry 



J K Block, N T Zinner and G M Bruun 

Department of Physics and Astronomy, University of Aarhus, DK-8000 Aarhus 
C, Denmark 

&mail: jkblock9phys.au.dk 

Abstract. We consider the density wave instability of fermionic dipoles aligned 

by an external field, and moving in equidistant layers at zero temperature. Using 
a conserving Hartree-Fock approximation, we show that correlations between 
dipoles in different layers significantly decrease the critical coupling strength for 
the formation of density waves when the distance between the layers is comparable 
to the inter-particle distance within each layer. This effect, which is strongest 
when the dipoles are oriented perpendicular to the planes, causes the density 
waves in neighboring layers to be in-phasc for all orientations of the dipoles. 
We furthermore demonstrate that the effects of the interlayer interaction can 
be understood from a classical model. Finally, we show that the interlayer 
correlations are important for experimentally relevant dipolar molecules, including 
the chemically stable ^'^Na.'^°K and '^°K^^^Cs, where the density wave regime is 
within experimental reach. 



PACS numbers: 03.75.Ss,67.85.d,68.65.Ac,73.20.Mf 
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1. Introduction 



The trapping and cooling of molecules in their rotational and vibrational ground state 
is a new research direction within the field of ultracold atomic and molecular physics 
[UminilllElinilZllHlini Uni • in contrast to the short-range isotropic interactions in 
typical cold atomic gases, molecules provide anisotropic potentials that typically have 
a long range dipolar part. This opens up a host of possibilities for exploring interesting 
physics [m [121 US] and also chemical reaction dynamics at low temperatures [El [15] . 

Chemical reaction losses can be large in three-dimensional (3D) samples [B] 
due to the attractive head-to-tail interaction of dipolar molecules. However, 
recent experiments using optical lattices [HI [IHj have shown that a low-dimensional 
confinement of the system can effectively suppress the loss. Interesting many-body 
phases have been proposed in such settings, including s- and p-wave superfluid states in 
single- and bilayer setups [Ml [171 UHl El EHl jST, '2T, '23 , density- waves in homogeneous 
[Ml US [Ml mi [2H1 HH] and lattice systems 130, iSl* i32j, and non-trivial Fermi liquid 
behavior [331 IMl ISS]- The long-range dipolar forces also opens up for a very rich 
spectrum of few-body bound state physics in ID [33 [33 133 [33] and 2D setups 
[10l[4Tl[42l[431|44]. 
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Figure 1. (a) Illustration of the setup considered. The dipoles reside in layers in 
the xj;-plane separated by the distance d. The angle 0, between the dipoles and 
the normal to the layers is in the xz-plane. The interaction between dipoles in the 
same layer, in adjacent layers, and in layers separated by a distance 2d is indicated 
by Vbi ^li and ^2 respectively, (b) The phase boundary between the normal (left) 
and the striped phase (right) in the (g, 6) plane for strictly 2D layers with ui = 0: 

a single layer ( ), two layers separated by dhPp = 0.5 ( ) and dfc^ = 1 

( — • — ), and three layers separated by dfc^ = 0.5 ( ) and dk^ = 1 ( — ■ — ). 

The ■ • line gives the RPA (Random Phase Approximation) result for a single 
layer for comparison. The case of a single quasi-2D layer with uifc^ = 0.1 is 
plotted by . 



We consider the density wave (stripe) instabilities of fermionic dipoles at zero 
temperature in 2D layers. The dipole moments are aligned by an external field, and 
they are moving in equidistant layers as illustrated in figure [ija). Calculating the 
density-density response function within a conserving Hartree-Fock approximation, 
we find the following. (I) The presence of several layers can decrease the critical 
coupling strength for stripe formation significantly. The effect is strongest for the 
dipoles oriented perpendicular to the planes whereas it vanishes when the dipoles 
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are aligned in the planes, see figure [T[b). (II) The mechanism for this decrease 
is the formation of stripes on top of each other in adjacent layers. This in-phase 
stripe pattern is always energetically favorable, independent of the orientation of the 
dipole moment which is somewhat surprising. It can be understood from a purely 
classical calculation which also explains the angle dependence of the effect. (Ill) 
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Figure 2. For a typical layer separation d = 1064 nm/2 and perpendicular 
polarization 8 = 0, the critical value of dipole moment times the square root of the 
mass, P\/rn, as a function of the square root of the two dimensional density for the 
single-, bi- and trilayer geometries. The horizontal lines indicate the permanent 
electrical dipole moment of the ground state (values from 45^) times the square 
root of the mass for the respective molecules. Outside the range of the axis lie 
^U^'^Rb and ^Li^^aQg ^j^j^ 40 g D^u and 64.6 Dy^ respectively. 

The decrease in the critical coupling strength for stripe formation due to interlayer 
correlations is significant for experimentally relevant systems consisting of ^Li^^K, 
^^Na^^K, 40k133Cs, '^U^^Rh and ^Li^SCs molecules. It is less important for ^^K^^Rb 
and ^Li^'^Na molecules which have smaller dipole moments. 

2. System 

The system consists of identical fermionic dipoles of mass m moving in equidistant 
layers separated by a distance d. Their dipole moment p is aligned forming the 
angle 9 with respect to the normal of the layers (z-axis) with their projection onto 
the planes defining the x-axis, see figure [T] The layers are formed by a deep ID 
optical lattice so that the dipoles in layer I {I an integer) reside in the lowest 
Wannier state in the z-direction. This is to a good approximation a Gaussian 
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fi{z) — exp[— (z — Id)"^ /2w^]7r^^/'^w^^/^ , where w is the width of the layer which 
is centered at Id. We neglect any trapping potential in the xy-plane so that the 
transverse states are labelled by the momentum k = {kx,ky) (we use units where 
h = kB = 1). 

In this basis the grand canonical Hamiltonian reads 
f k"^ \ 1 

^ ^^\^~ n^h^^^l^ 2^ ^ l^.i'(q)4+q,i4'-q,i'Ck',/'Ck,i, (1) 

k,i ^ ^ k,k',q 

where we have absorbed the harmonic oscillator energy of the z-direction in the 
chemical potential /z. Here, Ck_/ removes a dipole in layer I with momentum k. The 
interaction between two dipoles separated by r is 

where 9^ is the angle between r and the dipole moment p, and — /Attsq for electric 
dipoles. The effective interaction V;_//(q) between dipoles in layer I and I' is obtained 
by integrating the interaction V3d(?^ over the Gaussians ipi[z)ipii{z) combined with a 
2D Fourier transform. This yields (46j 



yo(q) = 



P2(cos0)-2^(^,^)F(q) 



(3) 



_3t(;v27r 

for the intralayer interaction where P2{x) = {Zx^ — l)/2 is the second Legendre 
polynomial, while 

F{q) = gexp[((jw)V2]erfc[(7w/%/2] and ^{9, ip) = cos{e)^ - sm{e)^ cos(cp)^ (4) 

The constant term in equation ([3| corresponds to a S{r) interaction which plays no 
role since we consider identical fermions. The real part of the interlayer interaction is 
only dependent on the difference / in layer numbers and is given by |47l I48j 

Vi{q) ^-27r D'a9,ip)qe'^''^ (5) 

for w <^ d. This approximation deviates less than 10% from the exact expression for 
w < d/5. In equation ([5]), we have only given the real part of the potential since the 
imaginary part is zero for momenta along the y-direction which are the ones relevant 
for stripe formation, as discussed in section [Sj 

The strength of the interaction is parametrized by the dimensionless number 

4mD^k°p . , 

where kp — \/Awri2D is defined from the 2D density. Likewise, a dimensionless measure 
for the layer separation is given by dk^p. The ratio of the layer distance to the typical 
inter-particle distance in a layer is dk'^p / y/in. 



3. Linear response and the Hartree-Fock approximation 

To analyze the instabilities of the homogeneous phase towards the formation of 
stripes, we consider the retarded density-density response function xf^ (r — r',t — t') = 
—i9{t — t'){[pi{r,t), pj{r' ,t')^) where denotes the layers and r = {x,y). The 

density operator for layer i is pi{r) — '0|(r)-0i(r) with ipi{r) the field operator for 
the dipoles in layer i. An instability toward the formation of a density wave with 
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wave number q shows up as a zero frequency = 0) pole of the Fourier transformed 
response function x^(q, t^)- 

We calculate the retarded response function using diagrammatic perturbation 
theory in Matsubara space Xii(ctc,*w„) with a;„ = 2mTT {n integer) a bosonic 
Matsubara frequency. The retarded function Xi5(qc,w) is then obtained by analytical 
continuation ia;„ — > w + iO"*" in the usual way |49j . As illustrated in figure |3j Xij (?) 
with q = (q, iw„) can be written as 

X^jlg) = 1] ['5u4,fc'n,(/c,Q) + n,(fc,g)r,,(fc,fc',g)n,(/c',g)] (7) 

k,k' 

where Ili{k,q) = Gi{k + q)Gi{k) is the particle-hole propagator with Gi{k) the single 
particle Green's function for the dipoles. The scattering matrix Tij{k,k' ,q), which 
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Figure 3. Top: Density-density response function. The thick lines indicate 
interacting Green's functions. Bottom: The Bethe-Salpeter equation for the 
particle- hole scattering matrix F. 



describes a particle with momentum k + q scattering on a hole with momentum k to 
produce a particle with momentum k' + q and a hole with momentum fc', obeys the 
Bethe-Salpeter equation 

r,y- (fc, k',q)^ ik,k',q)+J2 hi (fc, fc", q)Tll {k" , q)Tij (fc", k' , q) (8) 

l,k" 

as shown in figure [sj Here Iii{k,k" ,q) includes all scattering processes which are 
irreducible with respect to the particle-hole propagator 11; (A;, q). 

To proceed, we apply the self-consistent Hartree-Fock approximation illustrated 
in figure [4] All Green's functions are interacting in this approximation and the vertex 
Iii{k, k" , q) is given by the lowest order direct and exchange interactions. Writing the 
Green's function as G~[^{k) = ikn — £;k with ikn = (2n + 1)'kT a fermionic Matsubara 
frequency and ejk = k^/2TO — — we have the usual Hartree-Fock expression 

for the self-energy 

SKk) - ^[^/,j(0) - 5i,,Fo(k - k')]//k'. (9) 

Here /;k = [exp(/3e/k) + 1]^^ is the occupation of the k momentum state in layer I. 
In the calculations, we take /? = l/ksT — >■ oo appropriate for the zero temperature 
case. It follows from ([s]) that the Hartree shift Vij{0) for the energy of a dipole 
in layer I due to the interaction with the dipoles in layer j ^ I vanishes for thin 
layers, i.e. for w d. This is a consequence of the fact that the dipole-dipole 
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Figure 4. Hartree-Fock approximation for the Green's function and the 
irreducible particle-hole interaction. The thin lines indicate the non-interacting 
Green's functions and the wavy lines the interaction Vij . 



interaction integrates to zero over a plane as shown in [Appendix C[ In this paper, 
we keep w <C and since the density is the same in all layers, the self-energy and the 
chemical potential are both layer independent. We calculate the single particle Green's 
functions numerically, obtaining self-consistency through an iterative procedure. The 
scattering matrix Tij{k,k',q) is then determined from equation (Is]) using 



/,,,(k, k', q) = y,,(q) - 6,^,Vu{k - k') 



(10) 



Note that exchange is only included for interactions within the same layer as dipoles in 
different layers are distinguishable. Diagrammatically, this approximation corresponds 
to the summation of bubbles containing intralayer ladder interactions, which are 
connected to each other by inter- and intralayer interactions. 

Since the irreducible interaction vertex /^j is independent of frequency in this 
approximation, the particle-hole scattering matrix is independent of the internal 
frequencies ik'^,ik'^ and we can perform these frequency summations in equation ([s]), 
which then simplifies to 

ET ri T " \ /'k" — flk" 
luiK^ ,q) 



r„(k,k',g) -/„(k,k',q) 



iq„ 



-^r,,.(k",k',<z). 

£;k"+q 



(11) 



Equation (111 corresponds to a series of inhomogeneous Fredholm equations of the 
second kind in the first variable. From equation ([7| it follows that apart from 
the poles of 11 describing the particle-hole continuum, the density-density response 
function x and the particle-hole scattering matrix F share the same poles describing 
collective modes. Thus, we determine the critical coupling strength by searching for 
a zero frequency pole at a non-zero momentum q of the matrix TijCk, k', g), which is 
analytically continued to real frequency by ig„ to + iO^ . 

When 9 > 0, the anisotropic interaction favors dipoles that are aligned along 
the X-axis. This corresponds to a density wave with a wave vector pointing in the 
perpendicular direction ipc = 7r/2. Since the density wave is formed by particle-hole 
excitations, we expect the length of the wave vector to be ikpi^c) to minimize the 
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kinetic energy cost. This indeed follows directly from RPA calculations of the density- 
density response function [MJ I^S], and we expect it to hold even when exchange 
correlations are included. 

The self-consistent Hartree-Fock calculation of the response function is demanding 
numerically, and we describe in [Appendix A how it is implemented numerically. The 



payoff is that the approximation is conserving in the sense of Kadanoff-Baym [50] . We 
shall furthermore see that the exchange correlations have large effects on the critical 
coupling strength for the formation of stripes. 

4. Numerical results 

We now present numerical results for the critical coupling strength Qc as a function 
of dipole angle 9, layer separation dkp, and layer thickness wkp for fixed momentum 
in the direction ipc — with magnitude q = 2kp{ipc). We shall for concreteness 
consider the cases of one, two, and three layers. 

4-1. Single layer 

We first focus on the case of a single layer with vanishing thickness, i.e. w = Q. In 
figure [l|b) , we plot as a solid line the phase boundary between the normal phase 
(left) and the striped phase (right) for a single layer. The boundary has an intriguing 
non- monotonous behavior with a maximum critical coupling strength gc for 9 ~ 7r/4. 
For comparison we also plot as a dotted line the result of a RPA calculation using the 
interacting Green's functions [24 l I25 j . 

For small dipole angles 9, the RPA result underestimates the critical coupling 
strength significantly as compared to the conserving HF approximation. This 
demonstrates that exchange correlations suppress the formation of stripes. For larger 
angles, the shape of the phase boundary differs qualitatively from that obtained from 
the RPA calculation which predicts cx cos(6')^ (neglecting the effects of Fermi surface 
deformation). 

The dependency on 9 for the RPA result (green, dotted line of figure [ijb)) can be 
understood purely from the fact that the repulsive part of the interaction decreases as 
cos^ 9 as the dipoles are tilted towards the layer. For small 6', the exchange correlations 
suppress stripe formation. As 9 crosses the "magic angle" cos"^ {\ / \/?>) , the spatial (or 
momentum) average of the interaction goes from being repulsive to being attractive. 
Since exchange correlations enter through an average over the momentum transfer in 



the term ^(q) — V(\i — k'), sec ( 10 1 and the fc'-sum in ([8]), this means that the effects 
of the exchange term vanishes right at the magic angle as can be seen from figure [ijb). 
For larger 9, exchange correlations enhance the stripe instability and for 9 ^ tt/2 the 
instability is entirely driven by the exchange term since the direct interaction vanishes. 
This is the qualitative origin of the maximum in the critical coupling strength for an 
angle close to cos~^(l/\/3). 

Our result for the phase boundary agrees within 7% with that of ref. [55], and 
the critical coupling strength for = agrees within 5% of that reported in ref. |27) . 



As explained in [Appendix A we have taken care to use many /c-space points in the 



integration around the singular points in (11), and we estimate our results to be 
numerically accurate within 1%. On the other hand, our results differ substantially 
from those obtained using a self-consistently determined local field factor [29, . 
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We also plot in figure [IJb) the phase boundary for a quasi-2Z? single layer system 
where the finite depth of the ID optical lattice gives a width of wkp =0.1 for 
the Gaussian transverse wavef unctions. This softening of the z-direction reduces the 
repulsive part of the effective dipole-dipole interaction ([s]), and as a result the critical 
coupling strength for stripe formation is increased by up to 18%. For larger wkp, one 
has to take into account higher states in the ^-direction as investigated for the case of 
61 = in Ref. [22]. 

Note that we for simplicity have not included the region of p-wave superfluidity 
for 9 > 7r/4 [T^] in the phase diagram, nor the region of collapse due to a negative 
compressibility for large angles and coupling strengths [HI US [29j . 

4-2. Several layers 

In figure [ijb), we also plot the phase boundary in the case of two and three layers 
separated by the distances dkp — 0.5 and dkp — 1. This illustrates a main result of 
this paper: the presence of neighbouring layers reduces the critical coupling strength 
for stripe formation significantly when the layer distance d is comparable to the 
distance between particles within each layers. As expected, the effect decreases with 
increasing layer separation as follows directly from the exponential decay of the inter- 
layer interaction ([5|. This is illustrated further in figure [sja) where the critical 
coupling strength for the bi-layer case is compared to that of the single layer case 
as a function of layer separation and dipole angle. For small layer separation, the 
critical coupling strength can in fact be shown to scale as 1/N for ^ 1 layers, since 
the exchange correlations within each layer can be neglected in this limit, so that 
the problem reduces to that of dipoles with N internal degrees of freedom moving 
in a single layer [5B]. From figure [sjj a) we furthermore conclude that the effects of 
neighboring layers is strongest for small dipole angles. This effect has a simple classical 
interpretation as we shall demonstrate below. 

Figure [sjjb) shows the critical coupling strength for the bi-layer case for u; = 0. 
It increases with increasing layer separation dkp whereas it has an interesting angular 
dependence as a result of the interplay between the angular dependence of the 
interlayer and intralayer interaction. Note that for well-separated layers with w d, 




Figure 5. The critical coupling strengtli g^' for stripe formation for two layers 
as a function of dipole angle 9 and distance dk^ between the layers for w = 0: 
(a) gbi/g- and (b) g^i. 
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a non-zero layer thickness only changes the intra-layer interaction, and it therefore 
affects the critical coupling strength in a way similar to that of a single layer. 

It should be noted that s-wave interlayer superfluidity is likely to occur in the 
multilayer system. This has been discussed for perpendicular dipoles in for both 
bilayer [121 [23 HI] and multilayer cases (restricted to nearest-neighbour pairing) . 

For = 0, the interlayer superfluidity is expected for any value of g when density- 
waves are ignored. At sufficiently large g, we will likely have a competition of the two 
phases and a model containing both is necessary to infer which phase dominates or 
whether there can be coexistence. The superfluid could have p-wave symmetry for 
larger angles which may be more favorable to coexistence. This will be explored in 
future studies. 



4-3. Correlations between stripes in different layers 



To examine further the effects of neighboring layers, we now analyze the w = 
collective mode and the associated density oscillations at the critical coupling strength 
gc- We first focus on the case of two layers. The density oscillations induced by an 
external perturbation J^i J d'^rVi°^^{r, t)pi{r) are within linear response given as 



Spi{(i,uj) 




_Sp2{ci,Uj) 





Xia(q,^) 
X2a(q,w) 



xi.2(q,w) 
X2.2(q,w) 



(q>w) 



(12) 



At the critical coupling strength g^ one eigenvalue of the density response matrix in 



equation ( 12 1 diverges. We find that the mode which first diverges always is symmetric 
in the layer index I, independent of the dipole angle, except for 9 = Tr/2 where the 
modes are degenerate. Close to the critical coupling 1 — g/gc <^ 1, the density-density 
response matrix has the pole structure 



Xi,i 

X2,l 



Xl,2 
X2,2 



X 







x° 



1 - ff/ffc 



x'- 
x^ 



(13) 



as we demonstrate in detail in Appendix B Here, x is the Lindliard function including 



Hartree-Fock shifts of the single particle energies. Thus, the density instability in the 
two layers is in-phase and the stripes will be on top of each other as illustrated in 
figure ([T]). The mode which is anti-symmetric in the layer index becomes unstable at 
a larger coupling strength. This can be understood as a splitting of the eigen-mode 
for a single layer into a symmetric and anti-symmetric mode. 

The same result holds for more than two layers: The mode with the stripes in 
phase always becomes unstable first, irrespective of the value of 9, except for 9 = 7t/2 
where the modes are degenerate. This is illustrated in figure [6]ja) which depicts the 
critical coupling strength for the even and odd modes for the cases of two and three 
layers. The modes with higher g^ were calculated using the self-energy at the value for 
the lowest modes for simplicity. This changes the obtained values slightly but retains 
the relative ordering of the modes. For the case of three layers, there is an additional 
mode with no density fluctuations in the middle layer. It goes unstable for a coupling 
strength in between the values for the even and odd modes. This agrees with the 
results found within the RPA [53] [|} 



f The odd mode and the mode with no density modulations in the middle layer was mistakenly 
swapped around in the text of ref. I26| . 
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Figure 6. (a) The lowest density waves modes for the layer separation dk^ = 1 

and thickness w = 0. Bilayer in-phase ( ) and anti-phase ( ), and trilayer 

all layers in-phase ( ), the two outer layers in phase and the middle layer in 

anti-phase ( — • — ), and finally no density modulation in the middle layer and 

the outer layers in anti-phase ( ) (b) Classical interaction energy between two 

layers for the in-phase ( ) and anti-phase ( ) configurations as a function 

of the wavelength of the density modulations. 



5. Classical model 



We now demonstrate that the effects of neighboring layers can be understood from 
simple classical considerations. First, it is easy to show, see (C.2), that the classical 
interaction energy of a single dipole with an infinite layer of dipoles with homogenous 
density is zero. Second, we analyze the interaction between stripes by calculating the 
classical interaction energy between a single dipole and stripes of increased/decreased 
density of dipoles in a layer separated by a distance d. As explained in detail in 
[Appendix C[ a straightforward calculation gives that this interaction energy is 

Vclassical = T ^ ^l^^ ^ ^^^^^^ i^d/>^c) + SCCh^ (tTC?/ Ac)] , (14) 

where the — is for the in-phase case where one stripe in the plane is directly above 
the single dipole, and + is for the anti-phase case where the projection of the dipole 
onto the plane is in between two stripes. The linear density of dipoles within a stripe 
is given by 7, and Ac is the distance between the stripes within the layer. Equation 
( [l4| is plotted in figure |6][b). 

This demonstrates that the in-phase/out-of-phase configuration of the stripes has 
a negative/positive interaction energy thereby explaining why they decrease/increase 
the critical coupling strength Qc for stripe formation as compared to the single layer 
case. The cos^ 6* dependence in (14) furthermore explains why the effect decreases 
for increasing angle, with the modes becoming degenerate for = tt/2 as seen in 
figure [6]j a). It is reassuring that our full quantum mechanical calculation which is 
rather numerically involved, agrees with this classical analysis. 



6. Experimental realisations 

We now examine the importance of the interlayer correlations discussed above for 
typical experiments. In an experiment where the planar confinement is caused by a 
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ID optical lattice formed by counter propagation lasers with wavelength A, the layer 
distance will be d = A/2. Figure [2] shows the critical dipole strength times the square 
root of the mass, ^/mp, as a function of density for one, two, and three layers and 
the dipoles perpendicular to the planes for a typical wavelength of A = 1064 nm. 
For comparison, in [5] the JILA group reports a peak density y/ri2B = 0.58 /xm~^ or 
dkp = 1.1 in an experiment with ""^K^^Rb molecules in the rotational and vibrational 
ground state. In the figure, the horisontal lines are the permanent electrical dipole 
moment times the square root of the mass for several experimentally relevant fermionic 
dipolar molecules. Note that the effective dipole moment in the trap is somewhat 
smaller, since the dipoles are not aligned perfectly. It is however of the same order of 
magnitude as the permanent dipole moment; in [5] the JILA group reports a maximum 
value for the average dipole moment in the laboratory frame of about 40% of the 
permanent moment. 

From the figure, we see that the molecules ''Li^°K, ^^Na^'^K, and '*°K^'^'^Cs and 
moreover •^Li^^Rb and ^Li^^aQg^ which lie outside the fi gure, have such large values 
of y/mp, that one will observe stripe formation already in the regime of relatively low 
density where multilayer effects are important. This demonstrates the experimental 
relevance of the results discussed in this paper. The '*°K^''Rb and ^Li^'^Na molecules 
on the other hand have small values of \Jmp, and the density required to observe 
stripe formation is so high that interlayer correlations are not important. 

The molecules containing Li are all chemically unstable [5T] in the sense that 
the reaction 2YX — > Y2 + X2 is exothermic for Y=Li and X any other alkali metal. 
Molecules of ^'^Na^°K and ^°K^'^''Cs are however chemically stable in this sense, so 
they are prime candidates for studying the density wave instability. 

6.1. Experimental issues and detection 

Experiments with polar molecules operate at finite temperatures. The JILA 
experiments with KRb have reported temperatures down to T = 220 nK or T/Tp = 1.4 
[52j . This is close to but not quite in the degenerate regime, so a decrease in 
temperature is most likely needed to see many of the predicted phases. Furthermore, 
in the low-dimensional setups created by optical lattices, heating can occur as the 
lattice is turned on, demanding that the initial temperature be even lower to reach 
critical temperatures. 

In the strict 2D limit, the molecules are completely confined in the transverse 
direction. This means that their motion is reduced to a strict planar geometry. In 
addition, we do not allow any tunneling between the layers. The layer index can 
be considered an effective spin label. In this case, the non-zero temperature phase 
transition to ordered states such as the density wave is governed by the Berezinskii- 
Kosterlitz-Thouless (BKT) transition |53l|54j, and no true long-range order occurs. 
One consequence is that the BKT transition temperature is below the transition 
temperature obtained from mean-field theory. 

However, studies of one-dimensional arrays of tubes with two-component fermions 
that can undergo a pairing transition indicate that weak tunneling between the tubes 
can stabilize long-range order [55j. We imagine that this could also work for our 
multilayer system, i.e. by allowing weak tunneling between the layers the density- 
wave ordering becomes long-range and the transition temperature may be increased 
from the BKT value to the mean-field prediction. This intriguing possibility will be 
explored in future work. 
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For detection of the density wave state in the multilayered setup a number 
of different tecfiniques are possible. A quantum non-demolition measurement 
[551 1551 [571 1551 [5^ [BUI [5T] could be used to detect the large density fluctuations 
close to the phase transition. Alternatively, light scattering experiments proposed to 
detect dimerized pairing in multilayers |62j could be adapted to the density wave case. 

7. Conclusions 

Using a conserving Hartree-Fock approximation, we examined the density wave 
instability of aligned fermionic dipoles moving in equidistant planes. We found that 
while exchange correlations suppress the instability, it can be significantly enhanced 
by correlations between the layers. The inter-layer correlations exhibit an interesting 
dependence on the dipolar angle, and they result in the density waves in the different 
layers to be in-phase for all angles. We furthermore demonstrated that the physics 
of the interlayer correlations can be understood from a classical model. The density 
wave instability was shown to be experimentally accesible with realistic densities for 
experiments using ^Li^OR, ^Li^^Rb, ^Li^^SCs, ^^Na^"K and ^^K^^'^Cs molecules. For 
these molecules, interlayer correlations were furthermore predicted to decrease the 
critical coupling strength for density waves significantly. 
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Appendix A. Numerics 

Appendix A. 1. Green's function 

Since we use the conserving approximation, the limits on the self-energy integral in 
^ is dependent on S itself. Since we only consider T = 0, the ground state of the 
system is described by a deformed Fermi sea. Inspired by [28) . we make a Fourier 
series for the function h{ip) = kp{(p)/kp — 1 = '^n cos{2n(p) keeping the first 

six terms, which gives the deformation relative to the non-interacting Fermi sea, and 
determine the coefficients for a given interaction strength by an iterative procedure. 
An alternative is to approximate the Fermi sea by an ellipse and do a variational 
calculation as derived in [63] and applied in [16l [25]. Both approaches have been 
implemented by the authors of this paper, and the results are very similar. Interested 
readers can find the details of the procedure in [^5] . 



Appendix A. 2. Particle-hole scattering matrix 



The particle-hole scattering matrix is determined by (111, where the second variable 
(k') is not integrated, so it just appears as a parameter. 

The Fermi functions in ( 11 ) restrict the two dimensional integration domain to the 
(deformed) Fermi sea and a Fermi sea displaced by q. In order to reduce the numerical 
complexity, we define the functions r^(k,k',g) = rij(k, k',g) ±Tij{—'k— q, k',(7). 
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The inversion symmetry of both inter- and intralayer potentials V^(k) = y(— k) gives 
/(— k) = /(k) and e(— k — q) = £(k-|-q) which allows for a shift of the displaced Fermi 
sea to give for the static (w — 0) scattering matrix 

r± (k, k', q) = /± (k, k', q) + ^ /± (k, k", q) r± (k", k', q) , (A.l) 

where -^^j (k, k', q) = /ij(k, k', q) ± k — q, k', q). The scattering matrix is then 
given by r(k, k',q) — |[r+(k, k',q) + r^(k, k',q)] and thus has poles if and only if 
at least one of the symmetrized F's has a pole (unless the poles of r+ and happen 
to cancel). 

The Fermi sea is deformed by the dipole-dipole interaction, so we write the sum 
as an integral in polar coordinates and for fixed ip" parametrize the norm integral by 
x" = k" /kp{(p"). Let k be unit vector in the direction of k, then the integral becomes 

J24iux,<i)-^ — 

^ £ik" — £ik"+q 

k%^l + h{'p")x" 

' ei{k°p^l + hi^")x"k") - ei{kOp^l + hi^")x"k" + 2k%^l + h{n/2)ci) 



abscissa k^ and weights Wa for points in the singleTermi sea and approximate the 



An approach to solving the integral equation (|A.l|) for F^ is to choose suitable 



integral by a sum. Then (11) becomes a matrix equation where (i,j)'th block of the 
matrices describes the interaction on layer i from layer j. We suppress the common 
g-dependency and introduce the weights as a diagonal matrix W. In the double indices 
of layer number (roman letter) and k-grid point (Greek letter) the equation reads 

[5^a,iy - itj^Wi^xQ F± ^^.^ = /± (A.2) 

f . f II k I 

where the diagonal matrix x has entries Xi ■y — -eik + ' 

The irreducible interaction is finite, so the scattering matrix diverges when the 
matrix in the brackets becomes singular. This happens when the matrix Iia,i-yWijXi 



has an eigenvalue of 1 (see also Appendix B I . For / the direct interaction cancels 



7 



and the layers decouple as there is no exchange between different layers. For the single 
layer we find that the matrix has only negative eigenvalues, so F~ never contributes 
to the divergence. 

For fixed tpc — 7r/2 and q ~ 2kp{ipc) we vary g until the first eigenvalue crosses 1. 
The difference of the single particle energies of the particle and hole in the denominator 
of x^ gives rise to an integrable singularity at the edge of the integration region at 
(f" = 37r/2,fc" — 2kp{ipc), so we partition the (p" interval and cast 60 points in the 
[-7r/2-)-A(/j,37r/2-A(^] and 10 points in [37r/2-A((5, 37r/2+Av5]. For the x = k" /kF{p>) 
variable we cast 10 points in [0, 1 — Aa:] and 30 points in [1 — Ax, 1]. We choose the 
abscissa and weights in the four intervals by a Gauss-Legendre quadrature rule. For 
At/? = 0.05 and Ax = 0.02, a tripling of the number of points in either interval gives a 
change of the critical coupling strength that is less than the absolute tolerance 2 • 10"'^. 
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The matrix depends on g explicitly from the potential in /, but also indirectly 
through the self energies and Fermi function. The iterative procedure is as follows: 
Fo r the current gu ess for gc, calculate the deformation f{ip) of the Fermi sea according 



Appendix A.l 



rescale the fc-points and calculate the matrix 



to 

largest eigenvalue Ai and set the new guess to g = X^^ 
when the absolute change in g is less than 2 • 10^"^. 



'-Wi^x{j. find the 



The iteration is terminated 



Appendix B. Scattering matrix near the critical point 



As g ^ gc for fixed q Qc the matrix in the square brackets in (A. 2 ) becomes singular 



and the response blows up. If we make a spectral decomposition of i^^hXi-yi th'^ 
eigenvector with eigenvalue 1 will dominate in the vicinity of the divergence. To 
see this (for simplicity write / 



and for the diagonal matrix Wi^xjj = 
consider a g < gc, so that ID has eigenvalues A; < IVi. We need to find the inverse 
(1 — ID)~^, so we make an eigenvalue decomposition for ID, say ID — VAV^^ with 
A = diag(Ai, . . .). Then (ID)" = VA''V~^ ^ for n ^ oo. Now take the geometric 
series Ejr=o(-^^)" = ^EST^o^"^"^ = V"diag[l/(1 - Xi),...]V-^. So we have an 
inverse, as 



(1 - ID) y^ilDy = 1 - lim (/Z))" = 1. 



Tl = 



Inserting in (B.3 1 



r = l^diag[l/(l-Ai),...]F-i/ 
As (7 — ?> gc, the first 1/(1 — Ai) is divergent. 



(B.l) 



(B.2) 



Appendix B.l. Bilayer 

As mentioned above, the part is not divergent, so in the vicinity of the critical 
point, the scattering matrix is determined by the behaviour of r+, F k, . For 
the bilayer case we write the matrix structure of the 2n x 2n block matrix equation 



(A.2) explicitly as: the identity matrix Sa.^ = 1, the symmetric block matrices of 
the irreducible particle-hole interaction Iia,ij{ci = qc) = -^20, 27(^1 ^ ^c) = -^i", 
^ia,2'yi^ = qc) = l2a,i-yi^ = 1c) — I2 i diagonal product of the weights and 
particle-hole propagator WQ,x^(q = f\c)5a,-y = D and the particle-hole scattering 



matrix r+ j,g(q = qc) = r+ 



It 
It 



I^ 



D 




1.1 
2,1 



F 



2,2 



It 



^2 -'1 



(B.3) 



For g < gc we have the inversion in (B.2) above and now wish to examine the 



diagonalizing matrix V. The matrix I~^D has the symmetry found by switching blocks 

^01' 



by multiplying with C = 



1 



C has C ^ — C and commutes with both /+ and 



D since all 3 are (block) symmetric, thus it commutes with the product. C is it's 
own inverse, so it has eigenvalues ±1. Each has n linearly independent eigenvectors of 



the form v 



\w 



,±w ] for A = ±1 respectively. Because the matrices commute, the 



eigenvectors of I^D can be chosen to have the same form and the numerics show that 
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the eigenvector with the largest eigenvalue has the 
can be thus chosen to be of the block form 
'Wi W2 ' 

Wi -W2 



-form. The diagonalizing matrix 



V 



with V 



where the n x n matrices Wi are determined by the eigenvectors. As g 



(B.4) 
gc only the 

first 1/(1 — Ai) is divergent in the decomposition (B.2). To find the response functions, 
we need to sum all matrix elements in each block ([7|). Near the divergence all but the 
first entry in the diagonal matrix can be ignored, so 



2 



2(1-Ai) 



Vi 



(B.5) 



where vi 



Wi 



, and Wl is the first column of Wi from (|B.4|). From (|B.4|) we see, 

1 r 



that the first row of ^ has the structure — ^[zf , zf], where zf is the first row 
of Wl ^ . The matrix structure is now 



ri.i 
r2,i 



ri,2 
^2,2 



1 



wizlil+ + l+) w,zl{I+ + I+) 
wrzJ{I++I+) w^zf{I+ + I+) 



(B.I 



4(1 -Ai) 

We see that the blocks of the scattering matrix are all the same close to the divergence 
Tij = Fc. The r espo nse function is given by ([7|, so Xi,j is found by multiplying the 
block matrix in (B.6I by the diagonal matrix on both sides and then summing all 
matrix elements within each n x n block and finally adding the response function 
for non-interacting dipoles with HF single particle energy. The g-dependency of the 
/+I?-matrix comes both directly from the dipole-dipole interaction in /+ = gl^ where 
/+ is determined by the vectors k, k",q and from the self-energy in ■ 

The critical gc is defined by the nonlinear eigenvalue equation g{I^)D{g)v — Aiw, 
so by expanding near g = gc we see 1 — A oc — 5, so after the summation, the 2x2 
structure is 

(B.7) 



The numerics shows that the second largest eigenvalue value has an eigenvector 



Xi,i 


Xl,2 




■x° ■ 


1 




x"' 


.X2,l 


X2,2 




x°. 


1 - g/gc 


X 


x\ 



of the type V2 = [w^ 



Near the corresponding g, the inversion formula for 



{1 — (/+D)} does not work since {I'^D)'^ does not go to zero. This is because the 
system has already gone unstable. If we ignore this fact, say if the highest eigenvalue 
actually was the one with this eigenvector, the rank one approximation to the inverse 
would be the outer product of the (n -I- l)th column and row vectors of V and 

2 ^^2 : 



from (B.4 1 



T 



so the response near this gc2 would be 



Xi,i 


Xl,2 




\x' 1 


1 


■ xi 


-xi" 


.X2a 


X2,2 




x\ 


1 - g/gc 


.-xi 


X2 _ 



(B., 



Appendix C. Classical calculations 

First we show that a dipole in one layer does not feel a homogeneous distribution of 
dipoles in another (infinite) layer. The dipoles are in the xz-plane and aligned by the 
external field as p = D{sm-6E, 0, cos Oe) and the two layers are separated by a distance 
d. The angle between the relative vector and the dipole orientation is 

2a _ [P ■ r)"^ _ {xsinOE + dcos9E)^ 
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while dipole-dipole interaction is given by The calculations is done in a polar 
coordinate system in the xy-plane, — x'^ + , so 



I , I , , - (pcosipsineE + dcos0E 



= (C.2) 

This shows that any homogeneous background density cancels, so the interaction 
with a secondary layer is only dependent on the deviation from the average density. 
A density lower than the background can be modelled by changing the sign of the 
potential. 

Since the homogeneous background density does not contribute to the interlayer 
interaction, we can model the density wave phase by a periodic series of lines with 
alternating sign on the dipole-dipole interaction. The distance (wavelength of the 
density wave) between two lines with the same change in density is denoted Ac. The full 
quantum mechanical calculation ( 13 ) shows that the collective eigenmodes correspond 
to in-phase and anti-phase density modulations in the two layers, so we only consider 
these two extremes. 

We calculate the interaction between a single dipole in a stripe in one layer and 
the stripes in the other layer by splitting the lines into two sub-series. The first has 
lines directly on top while the lines of the other series are shifted by Ac/2. For the 
in-phase density modulations the first sub-series has excess density while the second 
has decreased density. For the anti-phase configuration the signs of the interactions 
are switched, but the geometry is otherwise the same. Using the coordinate system 
from figure [1} the stripes are parallel to the a;-axis. The relative coordinates between 
the dipole and a point in line n in the first sub-series and second sub-series is 

fi^n ^ [x,Xcn,d] f2,„ = [x, Ac(n- l/2),d], (C.3) 

respectively. For convenience we define — n^\^ + or a^^ ~ {n + ^)'^ + as 
the squared distance in the yz-plane in the two cases. The dipoles are assumed to 
be smeared out along the stripes with a linear dipole density of 7. So the interaction 
with line n is 

Summing the contributions from all stripes in both sub-series gives the interaction 
between a single dipole in a stripe with all the stripes in the other layer 

Vciassicai = T 1^^^^ [csch^(7rd/Ac) + sech^(7rd/Ac)] , (C.5) 

where the — is for the in-phase and + is for the anti-phase configuration. 
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